{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Nuclear Reactor \"Pincell\" Example\n",
    "In this example, we will demonstrate the creation of a simple 3D geometry. We will build a 3D \"pincell\" that is common in many pressurized water reactor geometries. It is basically a very long rectangular prism with a fuel cylinder runing through the middle of it lengthwise. The fuel cylinder is surrounded by a thin metal cladding region, with a very thin gap region between the fuel and the cladding. The area outside the cladding is the coolant. In this example, we will have the outer boundaries of the problems use a reflective boundary condition. The pincell will have the following dimensions:\n",
    "\n",
    "Height = 300 cm\n",
    "\n",
    "Pitch = 1.26 cm\n",
    "\n",
    "Fuel Radius = 0.4096 cm\n",
    "\n",
    "Inner Clad Radius = 0.4180 cm\n",
    "\n",
    "Outer Clad Radius = 0.4750 cm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Load the Library\n",
    "using ConstructiveSolidGeometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1 - Define Surfaces of the Problem\n",
    "\n",
    "Let's start with the rectangular prism boundaries (i.e., 6 planes).\n",
    "\n",
    "A **Plane** is defined by a point on the surface of the plane and the plane's unit normal vector. Both are specified in terms of a **Coord** object, which has three Float64 components x, y, and z corresponding to the vector's magnitude in each 3D cartesian direction.\n",
    "\n",
    "We can also select a boundary condition to assign to that surface that the user will be notified of at each collision. By default if no condition is specified, a transmission boundary is used allowing for a ray to pass through the surface without changing direction. Reflective or vacuum boundaries can also be specified by passing strings as the last argument of the surface constructor. The package will automatically enforce reflective boundary conditions. Vacuum boundaries behave identically to transmission boundaries, though the user may wish to enforce some other physics on the ray when a vacuum is encountered (or terminate the ray altogether)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "top =   Plane(Coord(0.0, 0.0, 150.0),  unitize(Coord(0.0, 0.0, 1.0)),  \"reflective\")\n",
    "bot =   Plane(Coord(0.0, 0.0, -150.0), unitize(Coord(0.0, 0.0, -1.0)), \"reflective\")\n",
    "left =  Plane(Coord(-.63, 0.0, 0.0),   unitize(Coord(-1.0, 0.0, 0.0)), \"reflective\")\n",
    "right = Plane(Coord(0.63, 0.0, 0.0),   unitize(Coord(1.0, 0.0, 0.0)),  \"reflective\")\n",
    "up =    Plane(Coord(0.0, 0.63, 0.0),   unitize(Coord(0.0, 1.0, 0.0)),  \"reflective\")\n",
    "down =  Plane(Coord(0.0, -0.63, 0.0),  unitize(Coord(0.0, -1.0, 0.0)), \"reflective\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's specify the infinite cylinder that will form the structure of the fuel, gap, and cladding regions.\n",
    "\n",
    "An **InfCylinder** is created by specifying a point along the central axis, the axial direction unit vector, and a radius. The point and directional vector are taken in the form of a **Coord** object. Note that these are transmission surfaces, as no boundary condition is specified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "clad_outer = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4750)\n",
    "clad_inner = InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4180)\n",
    "fuel =       InfCylinder(Coord(0.0, 0.0, 0.0), unitize(Coord(0.0, 0.0, 1.0)), 0.4096);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2  - Define the Halfspace Regions of the Problem\n",
    "Now we want to start building up our geometry using the surfaces we just defined. Each geometrical area, known as a **Cell**, is constructed logically out of any number of surfaces and their accompanying halfspaces. Surfaces divide all of space into two regions, positive and negative, according to the equation of the surface. For example, the negative halfspace of a sphere is the volume inside it, while the positive halfspace of a sphere is all space outside the sphere. We call each combination of a **Surface** and its halfspace sign a **Region**.\n",
    "\n",
    "We will begin by creating an empty array of **Cell** objects to store the new cells as we make them. This will eventually hold the water, cladding, gap, and fuel cells of our geometry. We also create an empty array of **Region** objects to hold the halfspace regions that define the water cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cells = Array{Cell}(0)\n",
    "regions = Array{Region}(0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by defining a set of regions which will be used to bound the water cell. A **Region** is defined as a surface and the sign of the halfspace that defines that region of space. I.e., a Region is created by a surface and a +1 or -1 integer. It is important to note that a plane may have two valid unit normals that are opposites of each other, but will affect the sign of the halfspace. Keep in mind that the (+) side of a plane's halfspace is on the side the normal points to."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(left, -1))\n",
    "push!(regions, Region(right, -1))\n",
    "push!(regions, Region(up, -1))\n",
    "push!(regions, Region(down, -1))\n",
    "push!(regions, Region(clad_outer, 1));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 - Define a Cell Using Logical Operators\n",
    "Now we define the logical manner in which we can combine the above regions to create the volume of the cell. In this library, we can use the intersection operator **^**, the union operator **|**, and the complement operator **~**. We define this logical operation using a Julia expression, which is wrapped by :() to indicate that it is an expression. Here, our halfspaces are defined so that we just want the simple intersection of all the regions. The numbers we used to identify the regions are the same as their indices in the regions array we just created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       ":(1 ^ (2 ^ (3 ^ (4 ^ (5 ^ (6 ^ 7))))))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ex = :(1 ^ 2 ^ 3 ^ 4 ^ 5 ^ 6 ^ 7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With that done, we can create the water cell and push it to our cells array that we made earlier. A new **Cell** is made by specifying the array of regions and the logical expression that defines the cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "push!(cells, Cell(regions, ex));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4 - Continue Making Cells\n",
    "Now we do the same thing for the gap, cladding, and water regions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Make the Cladding Cell\n",
    "regions = Array{Region}(0)\n",
    "\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(clad_outer, -1))\n",
    "push!(regions, Region(clad_inner, 1))\n",
    "\n",
    "ex = :(1 ^ 2 ^ 3 ^ 4)\n",
    "\n",
    "push!(cells, Cell(regions, ex))\n",
    "\n",
    "# Make the Gap Cell\n",
    "regions = Array{Region}(0)\n",
    "\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(clad_inner, -1))\n",
    "push!(regions, Region(fuel, 1))\n",
    "\n",
    "ex = :(1 ^ 2 ^ 3 ^ 4)\n",
    "\n",
    "push!(cells, Cell(regions, ex))\n",
    "\n",
    "# Make the Fuel Cell\n",
    "regions = Array{Region}(0)\n",
    "\n",
    "push!(regions, Region(top, -1))\n",
    "push!(regions, Region(bot, -1))\n",
    "push!(regions, Region(fuel, -1))\n",
    "\n",
    "ex = :(1 ^ 2 ^ 3)\n",
    "\n",
    "push!(cells, Cell(regions, ex));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5 - Store Cells In a Geometry Object\n",
    "A **Geometry** object holds an array of cells and a bounding **Box** that defines the scope of our problem. A **Box** object is an axis aligned 3D box defined by its lower left (minimum) corner and its upper right (maximum) corner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "bounding_box = Box(Coord(-.63, -.63, -150), Coord(.63, .63, 150))\n",
    "\n",
    "geometry = Geometry(cells, bounding_box);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6 - Plot a 2D Slice of the Geometry to Check Correctness\n",
    "To make sure that our geometry looks as intended, we will take a 2D slice of the problem space and plot it. This can be accomplished using the function **plot_geometry_2D(geometry::Geometry, view::Box, dim::Int64)** which takes a **Geometry** object, a view, and the pixel dimension resolution of the plot. The view is a 2D box (in x and y) that defines what will be plotted. The z dimension should be the same for lower_left and upper_right, and represents where the slice is taken."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_geometry_2D(geometry, Box(Coord(-0.63, -0.63, 0), Coord(.63, 0.63, 0)), 1000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also highlight a single cell within the geometry by using the function **plot_cell_2D(geometry::Geometry, view::Box, dim::Int64, cell_id::Int64)**, which is similar in function to the plot_geometry_2D() function, but also takes an additional argument indicating which cell to plot in black. Below, we plot the cladding cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src=\"\" />"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plot_cell_2D(geometry, Box(Coord(-0.63, -0.63, 0), Coord(.63, 0.63, 0)), 1000, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6 - Begin Ray Tracing\n",
    "At this point, we have a complete **Geometry** object that we can use to perform computations with. Let's test it by performing a standard ray trace operation. We utilize the generate_random_ray function to generate a **Ray** of random direction whose origin is uniformly distributed within a bounding box. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ray = generate_random_ray(geometry.bounding_box);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then find the point where the ray will intersect with next using the find_intersection function. This function returns a new **Ray** object that is at the location where the intersection occurs and with the updated direction if a reflection occured. The location has already been \"nudged\" slightly so that when the ray is next evaluated it will not result in an unpredictable evaluation due to the origin still being within machine precision distance of a surface. The function also returns the id of the surface that it intersected with as well as the boundary criteria of that surface. This is useful as in some cases one may wish to act upon the ray when it hits a reflective or vacuum boundary condition."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ray moved 0.06004697964739698 [cm] before hitting a transmission boundary\n"
     ]
    }
   ],
   "source": [
    "# Perform a single step of ray tracing on the geometry\n",
    "new_ray, id, boundary_type = find_intersection(ray, geometry)\n",
    "\n",
    "# Compute distance travelled by the ray\n",
    "distance = magnitude( new_ray.origin - ray.origin )\n",
    "\n",
    "println(\"Ray moved \", distance, \" [cm] before hitting a \", boundary_type, \" boundary\")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Julia 0.5.0",
   "language": "julia",
   "name": "julia-0.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
